%% fig4_DC_VI.m
clear; clc; close all;

% 忆阻器参数 (论文常取 a=0, b=1)
a = 0;
b = 1;

% 定义 phi 的扫描范围
phiRange = linspace(-3, 3, 1000);

% 根据 dphi/dt=0 => V = 0.2*phi^3 - phi
V = 0.2*phiRange.^3 - phiRange;

% 计算 I = (a + b*tanh(phi)) * V
I = (a + b*tanh(phiRange)) .* V;

% 计算 dI/dV (用于区分负斜率区)
% dI/dV = dI/dphi / dV/dphi
% 其中 dI/dphi = d/dphi[(a + b*tanh(phi))*V(phi)]
%             = ...
% dV/dphi = 0.2*(3*phi^2) - 1 = 0.6*phi^2 - 1
dVdphi = 0.6*phiRange.^2 - 1;
dWdphi = b * (1 - tanh(phiRange).^2); % 对于 (a + b*tanh(phi)) 的导数
dIdphi = dWdphi .* V + (a + b*tanh(phiRange)) .* (0.6*phiRange.^2 - 1) .* phiRange; 
% 上式分步写会更清晰
dIdV = dIdphi ./ dVdphi;  % 对应微分

% 将负斜率区与正斜率区分开
negIndex = find(dIdV < 0);
posIndex = find(dIdV >= 0);

% 绘图
figure; hold on;

% 绘制负斜率区(局部有源区)为红色
plot(V(negIndex), I(negIndex), 'r', 'LineWidth', 1.5);

% 绘制正斜率区为蓝色
plot(V(posIndex), I(posIndex), 'b', 'LineWidth', 1.5);

% 修饰
xlabel('V (V)');
ylabel('I (A)');
title('局部有源忆阻器的直流 V−I 特性图');
legend('负斜率区 (局部有源)','正斜率区','Location','best');
grid on;
axis equal;  % 使得坐标缩放一致
